% 在 6 组 (k, μ) 参数下绘制 MCLM 模型在 (x,y) 平面的相图。

clc; clear;

%参数设置
paramList = [  % 每行对应一组 (k, μ)
    0.03, 3.9; %   (a) (k, μ) = (0.03, 3.9)
    0.05, 3.6; %   (b) (k, μ) = (0.05, 3.6)
    0.25, 3.5; %   (c) (k, μ) = (0.25, 3.5)
    0.35, 3.6; %   (d) (k, μ) = (0.35, 3.6)
    0.35, 3.5; %   (e) (k, μ) = (0.35, 3.5)
    0.35, 3.1 %   (f) (k, μ) = (0.35, 3.1)
];

% 标签数组，用于标注子图 (a)~(f)
labels = {'(a)', '(b)', '(c)', '(d)', '(e)', '(f)'};

% 颜色数组（可根据论文调整颜色，这里示例用6种颜色）
colorList = { [0, 1, 0], [0.5, 0, 0.5], [1, 0.5, 0], [0, 0.5, 1], [1, 0, 1], [0.5, 0.5, 0.5] };

% % 初始条件
x0 = 0.1; 
y0 = 0; 
z0 = -1;

% 迭代设置
nTransient = 1000;   % 丢弃瞬态的步数
nPlot = 10000;        % 正式采样用于绘制相图的步数

for i = 1:6
    % 取当前参数 (k, μ)
    k = paramList(i, 1);
    mu = paramList(i, 2);
    
    % 初始化状态
    state = [x0, y0, z0];
    
    % (a) 丢弃瞬态迭代：使系统进入长期行为区域
    for n = 1:nTransient
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
    end
    
    % (b) 正式采样：采样 nPlot 步，将每步的 x 值记录下来
    xData = zeros(1, nPlot);
    yData = zeros(1, nPlot);
    for n = 1:nPlot
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
        xData(n) = dx;
        yData(n) = dy;
    end
    
    % 在 3×2 布局中，第 i 个子图
    subplot(3,2,i);
    scatter(xData, yData, 1, colorList{i}, 'filled'); 

    % 固定 x,y 范围为 [0,1.2]
    xlim([0,1.2]);
    ylim([0,1.2]);
    axis equal;  % 保持 x,y 比例相同
    grid on;

    xlabel('x'); 
    ylabel('y');

end
